############################################
## authors:   Michael L. Wicki & Robert A. Huber
## contact:   michael.wicki@istp.ethz.ch, ETH Zurich
## file name: JPP_replication_code.R
## Context:   ERC Mobility Project on Public Support and Policy Packaging
## started:   2018-07-19
## Summary:   Analyses the main results of the JPP Paper entitled "Can Policy-Packaging Increase Public Support for Costly Policies? Insights from a Choice Experiment on Policies against Vehicle Emissions"
##############################################

#set working directory if necessary
#setwd()

#load dataset
load("replication_data.RData")

#load necessary packages
library(cjoint)
library(data.table)
library(dplyr)
library(nFactors)
library(ggplot2)
library(lme4)
library(effects)
library(ggthemes)
library(tidyr)
library(forcats)
library(texreg)
library(influence.ME)


######################## Figure 2 ############################
#models

m_base_choice <- lmer(choice ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                        attrib5_lab +  attrib6_lab  +  attrib8_lab +  
                        attrib9_lab +  attrib7_lab + (1|id), data = df_conj_out2)

m_base_rate <- lmer(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                      attrib5_lab +  attrib6_lab +  attrib8_lab +  
                      attrib9_lab +  attrib7_lab + (1|id), data = df_conj_out2)

m_control_rate <- lmer(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                      attrib5_lab +  attrib6_lab +  attrib8_lab +  
                      attrib9_lab +  attrib7_lab + frame + age + gender + edu + usagecar + (1|id), data = df_conj_out2)


#library(base)
c <- c(as.numeric(coef(m_base_choice)$id[1,]), as.numeric(coef(m_base_rate)$id[1,]))
s <- c(as.numeric(se.fixef(m_base_choice)), as.numeric(se.fixef(m_base_rate)))

#remove control coef.
var <- factor(rep(c("Intercept", "Registration Ban", "Stricter Energy Label",
                     "Road Pricing",
                     "Emission Dependent Car Tax", "Road and PT Infrastructure",
                     "General Budget", 
                     "PT Price Reduction",
                     "Campaign", "2030", "2035", "2040", "2045",
                    "Test Phase", "Expert Involvement"), 2),
         levels = c("Intercept", "Registration Ban", "Stricter Energy Label",
                    "Road Pricing",
                    "Emission Dependent Car Tax", "Road and PT Infrastructure",
                    "General Budget", 
                    "PT Price Reduction",
                    "Test Phase", "2030", "2035", "2040", "2045", "Expert Involvement",
                    "Campaign"))

var <- factor(var, levels = c("Intercept", "Registration Ban", "Stricter Energy Label",
                              "Road Pricing",
                              "Emission Dependent Car Tax", "Road and PT Infrastructure",
                              "General Budget", 
                              "PT Price Reduction",
                              "Campaign",
                              "Test Phase", "Expert Involvement", "2030", "2035", "2040", "2045"))

var <- fct_rev(var)

dv <- factor(rep(c("Choice", "Rating"), each=15))

df_conj_plot <- data.frame(c,s,var,dv)

df_conj_plot$dv_dodge <- factor(ifelse(df_conj_plot$dv == "Choice", "Choice", "Rating"))

df_conj_plot$primary <- factor(ifelse(df_conj_plot$var == "Registration Ban" | df_conj_plot$var == "Stricter Energy Label" | df_conj_plot$var == "Road Pricing" | df_conj_plot$var == "Emission Dependent Car Tax", "Primary Measures",
                                      ifelse(df_conj_plot$var == "Campaign" | df_conj_plot$var == "Test Phase" | df_conj_plot$var == "Expert Involvement", "Ancillary Measures", "Other Policy Features")))

df_conj_plot$primary <- fct_rev(df_conj_plot$primary)

facet_grid(policy~grp, scales = "free")
  
df_conj_plot$policy <- ifelse(df_conj_plot$var == "Registration Ban", 1, ifelse(
  df_conj_plot$var == "Stricter Energy Label", 2, ifelse(
    df_conj_plot$var == "Road Pricing", 3, ifelse(
      df_conj_plot$var == "Emission Dependent Car Tax", 4, ifelse(
        df_conj_plot$var == "Road and PT Infrastructure" | df_conj_plot$var == "General Budget" | df_conj_plot$var == "PT Price Reduction", 5, ifelse(
          df_conj_plot$var == "Campaign", 6, ifelse(
            df_conj_plot$var == "Test Phase", 7, ifelse(
              df_conj_plot$var == "Expert Involvement", 8, ifelse(
                df_conj_plot$var == "2030" | df_conj_plot$var == "2035" | df_conj_plot$var == "2040" | df_conj_plot$var == "2045", 9, 
                NA)))))))))

df_conj_plot$policy = factor(df_conj_plot$policy, levels=c(1:9), labels=c("Primary\nMeasures","Primary\nMeasures","Primary\nMeasures", "Primary\nMeasures", "Revenue\nUse", "Ancillary\nMeasures", "Ancillary\nMeasures", "Ancillary\nMeasures", "Implementation\nYear"))


#Export Figure 2
p_base <- ggplot(subset(df_conj_plot, var != "Intercept"), aes(x=var, y=c, group = primary)) +
  geom_pointrange(aes(ymin=c-1.96*s, ymax= c+1.96*s, colour= primary, shape = primary), position=position_dodge(width=1), size=1)+
  theme_bw() + 
  coord_flip() +
  ylab("Change in E[Y]") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") + 
  facet_grid(policy~dv_dodge, scales = "free") +
  theme(strip.background = element_rect(fill = 'white'))+
  scale_colour_manual(name = "Measure Type", 
                      values = c("#d01c8b", "#b8e186", "#4dac26"),  
                      labels = c("Primary Measures", "Other Policy Features", "Ancillary Measures"))+
  scale_shape_manual(name = "Measure Type",
                     values = c(15, 17, 19),
                     labels = c("Primary Measures", "Other Policy Features", "Ancillary Measures"))+
  theme(text = element_text(size=16),
        legend.position="bottom")

ggsave(p_base, filename = "./figs/packaging/JPP/JPP_Figure2.png", width = 10.5, height = 7.5, unit="in", dpi = 600)

######################## Figure 4 ############################

df_conj_out2$sumAnc  <- ifelse(df_conj_out2$Attrib6 == "_2", 1, 0) +
  ifelse(df_conj_out2$Attrib8 == "_2", 1, 0) + 
  ifelse(df_conj_out2$Attrib9 == "_2", 1, 0)

m_choice_anc <- lmer(choice ~ attrib1_lab*sumAnc +  attrib2_lab*sumAnc +  attrib3_lab*sumAnc +  attrib4_lab*sumAnc +  
                        attrib5_lab +  attrib7_lab + I(sumAnc^2) +  (1|id), data = df_conj_out2)

m_rate_anc <- lmer(rate ~ attrib1_lab*sumAnc +  attrib2_lab*sumAnc +  attrib3_lab*sumAnc +  attrib4_lab*sumAnc +  
                       attrib5_lab +  attrib7_lab + I(sumAnc^2) +  (1|id), data = df_conj_out2)

m_rate_control_anc <- lmer(rate ~ attrib1_lab*sumAnc +  attrib2_lab*sumAnc +  attrib3_lab*sumAnc +  attrib4_lab*sumAnc +  
                     attrib5_lab +  attrib7_lab + I(sumAnc^2) + frame + age + gender + edu + usagecar + (1|id), data = df_conj_out2)


#Choice###

fit1 <- effect("attrib1_lab*sumAnc", m_choice_anc, xlevels = list(sumAnc=c(0:3)), 
              x.var="sumAnc", confidence.level = 0.95)
library(plyr)

fit1 <- data.frame(fit1)

fit1 <- rename(fit1, c("attrib1_lab"="prim"))


fit2 <- effect("sumAnc*attrib2_lab", m_choice_anc, xlevels = list(sumAnc=c(0:3)), 
               x.var="sumAnc", confidence.level = 0.95)

fit2 <- data.frame(fit2)

fit2 <- rename(fit2, c("attrib2_lab"="prim"))

fit3 <- effect("sumAnc*attrib3_lab", m_choice_anc, xlevels = list(sumAnc=c(0:3)), 
               x.var="sumAnc", confidence.level = 0.95)

fit3 <- data.frame(fit3)

fit3 <- rename(fit3, c("attrib3_lab"="prim"))

fit4 <- effect("sumAnc*attrib4_lab", m_choice_anc, xlevels = list(sumAnc=c(0:3)), 
               x.var="sumAnc", confidence.level = 0.95)

fit4 <- data.frame(fit4)

fit4 <- rename(fit4, c("attrib4_lab"="prim"))



df_sumAnc <- rbind(fit1, fit2, fit3, fit4)

df_sumAnc$subsetPrim <- factor(ifelse(grepl(x = df_sumAnc$prim, pattern = "1)"), "remove", "keep"))

df_sumAnc$sumAnc <- factor(df_sumAnc$sumAnc)

df_sumAnc_choice <- df_sumAnc

# Rate ###

fit1 <- effect("attrib1_lab*sumAnc", m_rate_anc, xlevels = list(sumAnc=c(0:3)), 
               x.var="sumAnc", confidence.level = 0.95)
library(plyr)

fit1 <- data.frame(fit1)

fit1 <- rename(fit1, c("attrib1_lab"="prim"))


fit2 <- effect("sumAnc*attrib2_lab", m_rate_anc, xlevels = list(sumAnc=c(0:3)), 
               x.var="sumAnc", confidence.level = 0.95)

fit2 <- data.frame(fit2)

fit2 <- rename(fit2, c("attrib2_lab"="prim"))

fit3 <- effect("sumAnc*attrib3_lab", m_rate_anc, xlevels = list(sumAnc=c(0:3)), 
               x.var="sumAnc", confidence.level = 0.95)

fit3 <- data.frame(fit3)

fit3 <- rename(fit3, c("attrib3_lab"="prim"))

fit4 <- effect("sumAnc*attrib4_lab", m_rate_anc, xlevels = list(sumAnc=c(0:3)), 
               x.var="sumAnc", confidence.level = 0.95)

fit4 <- data.frame(fit4)

fit4 <- rename(fit4, c("attrib4_lab"="prim"))



df_sumAnc <- rbind(fit1, fit2, fit3, fit4)

df_sumAnc$subsetPrim <- factor(ifelse(grepl(x = df_sumAnc$prim, pattern = "1)"), "remove", "keep"))

df_sumAnc$sumAnc <- factor(df_sumAnc$sumAnc)

df_sumAnc_rate <- df_sumAnc

df_sumAnc <- rbind(df_sumAnc_choice, df_sumAnc_rate)

df_sumAnc$dv <- factor(rep(c("Choice", "Rating"), each=32))

df_sumAnc$dv_dodge <- factor(ifelse(df_sumAnc$dv == "Choice", "Choice", "Rating"))

df_sumAnc$yint <- ifelse(df_sumAnc$dv == "Choice", 0.5, 4)

df_sumAnc$prim <- substring(df_sumAnc$prim, 4)

df_sumAnc$prim <- factor(df_sumAnc$prim, levels = c("Car Ban", "Tighter Energy Label", "Road Pricing", "Car Tax"),
                         labels = c(c("Registration\nBan", "Stricter\nEnergy\nLabel", "Road\nPricing", "Emission\nDependent\nCar Tax")))

df_sumAnc$ylabel <- ifelse(df_sumAnc$dv == "Choice", .65, 4.5)
df_sumAnc$ylabelMargins <- ifelse(df_sumAnc$dv == "Choice", 0.13, 0.29)

#Figure A3
p_sumAnc <- ggplot(subset(df_sumAnc, subsetPrim == "keep"), aes(x= prim, y=fit, group = sumAnc, shape = sumAnc))+
  geom_pointrange(aes(ymin=lower, ymax= upper), position=position_dodge(width=.5), size=1)+
  theme_bw() + 
  facet_wrap(~dv_dodge, scales = "free_y") +
  ylab("Predicted Support") +
  xlab("Primary Measures") +
  geom_hline(aes(yintercept = yint), lty="dashed") +
  scale_shape_discrete(name = "Number of\nAncillary\nMeasures") +
  theme(strip.background = element_rect(fill = 'white')) +
  theme(text = element_text(size=16)) +
  geom_vline(xintercept = 2.5, lty="dashed") +
  geom_text(aes(y = as.numeric(ylabel), x=2.4),label = ('Command-and-Control'), hjust = 1)+
  geom_text(aes(y = as.numeric(ylabel), x=2.6), label = ('Market-based'), hjust = 0)

#Margins ###

df_sumAnc$relFit <- ifelse(df_sumAnc$prim == "Registration\nBan" & df_sumAnc$dv == "Choice", df_sumAnc$fit - df_sumAnc$fit[which(df_sumAnc$prim == "Registration\nBan" & df_sumAnc$sumAnc == 0 & df_sumAnc$dv == "Choice")],
                     ifelse(df_sumAnc$prim == "Registration\nBan" & df_sumAnc$dv == "Rating", df_sumAnc$fit - df_sumAnc$fit[which(df_sumAnc$prim == "Registration\nBan" & df_sumAnc$sumAnc == 0 & df_sumAnc$dv == "Rating")],
                            ifelse(df_sumAnc$prim == "Emission\nDependent\nCar Tax" & df_sumAnc$dv == "Choice", df_sumAnc$fit - df_sumAnc$fit[which(df_sumAnc$prim == "Emission\nDependent\nCar Tax" & df_sumAnc$sumAnc == 0 & df_sumAnc$dv == "Choice")],
                                   ifelse(df_sumAnc$prim == "Emission\nDependent\nCar Tax" & df_sumAnc$dv == "Rating", df_sumAnc$fit - df_sumAnc$fit[which(df_sumAnc$prim == "Emission\nDependent\nCar Tax" & df_sumAnc$sumAnc == 0 & df_sumAnc$dv == "Rating")],
                                          ifelse(df_sumAnc$prim == "Stricter\nEnergy\nLabel" & df_sumAnc$dv == "Choice", df_sumAnc$fit - df_sumAnc$fit[which(df_sumAnc$prim == "Stricter\nEnergy\nLabel" & df_sumAnc$sumAnc == 0 & df_sumAnc$dv == "Choice")],
                                                 ifelse(df_sumAnc$prim == "Stricter\nEnergy\nLabel" & df_sumAnc$dv == "Rating", df_sumAnc$fit - df_sumAnc$fit[which(df_sumAnc$prim == "Stricter\nEnergy\nLabel" & df_sumAnc$sumAnc == 0 & df_sumAnc$dv == "Rating")],
                                                        ifelse(df_sumAnc$prim == "Road\nPricing" & df_sumAnc$dv == "Choice", df_sumAnc$fit - df_sumAnc$fit[which(df_sumAnc$prim == "Road\nPricing" & df_sumAnc$sumAnc == 0 & df_sumAnc$dv == "Choice")],
                                                               ifelse(df_sumAnc$prim == "Road\nPricing" & df_sumAnc$dv == "Rating", df_sumAnc$fit - df_sumAnc$fit[which(df_sumAnc$prim == "Road\nPricing" & df_sumAnc$sumAnc == 0 & df_sumAnc$dv == "Rating")], NA))))))))


p_sumAnc_margins <- ggplot(subset(df_sumAnc, subsetPrim == "keep"), aes(x= prim, y=relFit, group = sumAnc, shape = sumAnc))+
  geom_pointrange(aes(ymin=relFit - 1.96*se, ymax= relFit + 1.96*se), position=position_dodge(width=.5), size=1)+
  theme_bw() + 
  facet_grid(~dv_dodge, scales = "free_y") +
  ylab("Predicted Support") +
  xlab("Primary Measures") +
  geom_hline(yintercept = 0, lty="dashed") +
  scale_shape_discrete(name = "Number of\nAncillary\nMeasures") +
  theme(strip.background = element_rect(fill = 'white')) +
  theme(text = element_text(size=16))+
  geom_vline(xintercept = 2.5, lty="dashed")  +
  annotate("text", x = 2.4, y = .3, label = c('Command-and-Control'), hjust = 1)+
  annotate("text", x = 2.6, y = .3, label = c('Market-based'), hjust = 0)

###### Figure 4 export
ggsave(p_sumAnc_margins, filename = "./figs/packaging/JPP/JPP_Figure4.png", width = 12, height = 7, unit="in", dpi = 600)

######################## Figure 3 ############################
# Rating ###

mCB_rateI4 <- lmer(rate ~ attrib1_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

mEL_rateI4 <- lmer(rate ~ attrib2_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

mRP_rateI4 <- lmer(rate ~ attrib3_lab *attrib6_lab * attrib8_lab*  attrib9_lab  + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

mCT_rateI4 <- lmer(rate ~ attrib4_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

# Car Ban ###


fit <- effect("attrib1_lab:attrib6_lab:attrib8_lab:attrib9_lab", mCB_rateI4, xlevels = list(attrib1_lab=c("1) No Car Ban","2) Car Ban"),
                                                                                            attrib6_lab=c("1) No Campaign", "2) Campaign"),
                                                                                            attrib8_lab=c("1) No Test Phase", "2) Test Phase"),
                                                                                            attrib9_lab=c("1) Fedeal Gov't", "2) Federal Gov't and Experts")), 
              x.var="attrib1_lab", confidence.level = 0.95)

fit1 <- data.frame(fit)

fit1 <- subset(fit1, attrib1_lab == "2) Car Ban")

fit1$labels <- c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement")

fit1$labels <- factor(fit1$labels, levels = c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                                              "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                                              "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement"))

fit1 <- rename(fit1, c("attrib1_lab" = "prim"))

# Energy Label ###

fit <- effect("attrib2_lab:attrib6_lab:attrib8_lab:attrib9_lab", mEL_rateI4, xlevels = list(attrib2_lab=c("1) No Tighter Energy Label","2) Tighter Energy Label"),
                                                                                            attrib6_lab=c("1) No Campaign", "2) Campaign"),
                                                                                            attrib8_lab=c("1) No Test Phase", "2) Test Phase"),
                                                                                            attrib9_lab=c("1) Fedeal Gov't", "2) Federal Gov't and Experts")), 
              x.var="attrib2_lab", confidence.level = 0.95)

fit2 <- data.frame(fit)

fit2 <- subset(fit2, attrib2_lab == "2) Tighter Energy Label")

fit2$labels <- c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                 "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                 "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement")

fit2$labels <- factor(fit2$labels, levels = c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                                              "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                                              "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement"))

fit2 <- rename(fit2, c("attrib2_lab" = "prim"))

## Road Pricing ###

fit <- effect("attrib3_lab:attrib6_lab:attrib8_lab:attrib9_lab", mRP_rateI4, xlevels = list(attrib3_lab=c("1) No Road Pricing","2) Road Pricing"),
                                                                                            attrib6_lab=c("1) No Campaign", "2) Campaign"),
                                                                                            attrib8_lab=c("1) No Test Phase", "2) Test Phase"),
                                                                                            attrib9_lab=c("1) Fedeal Gov't", "2) Federal Gov't and Experts")), 
              x.var="attrib3_lab", confidence.level = 0.95)

fit3 <- data.frame(fit)

fit3 <- subset(fit3, attrib3_lab == "2) Road Pricing")

fit3$labels <- c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                 "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                 "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement")

fit3$labels <- factor(fit3$labels, levels = c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                                              "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                                              "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement"))

fit3 <- rename(fit3, c("attrib3_lab" = "prim"))

# Car Tax ###


fit <- effect("attrib4_lab:attrib6_lab:attrib8_lab:attrib9_lab", mCT_rateI4, xlevels = list(attrib4_lab=c("1) No Car Tax","2) Car Tax"),
                                                                                            attrib6_lab=c("1) No Campaign", "2) Campaign"),
                                                                                            attrib8_lab=c("1) No Test Phase", "2) Test Phase"),
                                                                                            attrib9_lab=c("1) Fedeal Gov't", "2) Federal Gov't and Experts")), 
              x.var="attrib4_lab", confidence.level = 0.95)

fit4 <- data.frame(fit)

fit4 <- subset(fit4, attrib4_lab == "2) Car Tax")

fit4$labels <- c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                 "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                 "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement")

fit4$labels <- factor(fit4$labels, levels = c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                                              "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                                              "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement"))

fit4 <- rename(fit4, c("attrib4_lab" = "prim"))

fit_rate <- rbind(fit1, fit2, fit3, fit4)

# Choice ###

mCB_choiceI4 <- lmer(choice ~ attrib1_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

mEL_choiceI4 <- lmer(choice ~ attrib2_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

mRP_choiceI4 <- lmer(choice ~ attrib3_lab *attrib6_lab * attrib8_lab*  attrib9_lab  + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

mCT_choiceI4 <- lmer(choice ~ attrib4_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

# Car Ban ###


fit <- effect("attrib1_lab:attrib6_lab:attrib8_lab:attrib9_lab", mCB_choiceI4, xlevels = list(attrib1_lab=c("1) No Car Ban","2) Car Ban"),
                                                                                            attrib6_lab=c("1) No Campaign", "2) Campaign"),
                                                                                            attrib8_lab=c("1) No Test Phase", "2) Test Phase"),
                                                                                            attrib9_lab=c("1) Fedeal Gov't", "2) Federal Gov't and Experts")), 
              x.var="attrib1_lab", confidence.level = 0.95)

fit1 <- data.frame(fit)

fit1 <- subset(fit1, attrib1_lab == "2) Car Ban")

fit1$labels <- c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                 "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                 "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement")

fit1$labels <- factor(fit1$labels, levels = c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                                              "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                                              "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement"))

fit1 <- rename(fit1, c("attrib1_lab" = "prim"))

# Energy Label ###

fit <- effect("attrib2_lab:attrib6_lab:attrib8_lab:attrib9_lab", mEL_choiceI4, xlevels = list(attrib2_lab=c("1) No Tighter Energy Label","2) Tighter Energy Label"),
                                                                                            attrib6_lab=c("1) No Campaign", "2) Campaign"),
                                                                                            attrib8_lab=c("1) No Test Phase", "2) Test Phase"),
                                                                                            attrib9_lab=c("1) Fedeal Gov't", "2) Federal Gov't and Experts")), 
              x.var="attrib2_lab", confidence.level = 0.95)

fit2 <- data.frame(fit)

fit2 <- subset(fit2, attrib2_lab == "2) Tighter Energy Label")

fit2$labels <- c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                 "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                 "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement")

fit2$labels <- factor(fit2$labels, levels = c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                                              "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                                              "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement"))

fit2 <- rename(fit2, c("attrib2_lab" = "prim"))

## Road Pricing ###

fit <- effect("attrib3_lab:attrib6_lab:attrib8_lab:attrib9_lab", mRP_choiceI4, xlevels = list(attrib3_lab=c("1) No Road Pricing","2) Road Pricing"),
                                                                                            attrib6_lab=c("1) No Campaign", "2) Campaign"),
                                                                                            attrib8_lab=c("1) No Test Phase", "2) Test Phase"),
                                                                                            attrib9_lab=c("1) Fedeal Gov't", "2) Federal Gov't and Experts")), 
              x.var="attrib3_lab", confidence.level = 0.95)

fit3 <- data.frame(fit)

fit3 <- subset(fit3, attrib3_lab == "2) Road Pricing")

fit3$labels <- c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                 "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                 "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement")

fit3$labels <- factor(fit3$labels, levels = c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                                              "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                                              "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement"))

fit3 <- rename(fit3, c("attrib3_lab" = "prim"))

# Car Tax ###


fit <- effect("attrib4_lab:attrib6_lab:attrib8_lab:attrib9_lab", mCT_choiceI4, xlevels = list(attrib4_lab=c("1) No Car Tax","2) Car Tax"),
                                                                                            attrib6_lab=c("1) No Campaign", "2) Campaign"),
                                                                                            attrib8_lab=c("1) No Test Phase", "2) Test Phase"),
                                                                                            attrib9_lab=c("1) Fedeal Gov't", "2) Federal Gov't and Experts")), 
              x.var="attrib4_lab", confidence.level = 0.95)

fit4 <- data.frame(fit)

fit4 <- subset(fit4, attrib4_lab == "2) Car Tax")

fit4$labels <- c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                 "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                 "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement")

fit4$labels <- factor(fit4$labels, levels = c("Base Line", "w/ Campaign", "w/ Test Phase", "w/ Campaign & Test Phase",
                                              "w/ Expert Involvement", "w/ Campaign & Expert Involvement",
                                              "w/ Test Phase & Expert Involvement", "w/ Campaign, Test Phase & Expert Involvement"))

fit4 <- rename(fit4, c("attrib4_lab" = "prim"))

fit_choice <- rbind(fit1, fit2, fit3, fit4)

fit <- rbind(fit_choice, fit_rate)

fit$dv <- factor(rep(c("Choice", "Rating"), each = 32))

fit$prim <- substring(fit$prim, 4)

fit$yint <- ifelse(fit$dv == "Choice", 0.5, 4)

fit$col <- factor(ifelse(grepl(x = fit$attrib6_lab, pattern = "1)"), 0, 1) +
  ifelse(grepl(x = fit$attrib8_lab, pattern = "1)"), 0, 1) +
  ifelse(grepl(x = fit$attrib9_lab, pattern = "1)"), 0, 1))
  
fit$prim <- factor(fit$prim, levels = c("Car Ban", "Tighter Energy Label", "Road Pricing", "Car Tax"),
                   labels = c(c("Registration\nBan", "Stricter\nEnergy\nLabel", "Road\nPricing", "Emission\nDependent\nCar Tax")))

fit$labels <- factor(fit$labels, levels = c("Base Line", "w/ Campaign", "w/ Expert Involvement", "w/ Test Phase",
                                            "w/ Campaign & Expert Involvement", "w/ Campaign & Test Phase", "w/ Test Phase & Expert Involvement",
                                            "w/ Campaign, Test Phase & Expert Involvement"))

p_Anc_choice <- ggplot(subset(fit, dv == "Choice"), aes(x= prim, y=fit, group = labels, shape = labels, colour = col, linetype = col))+
  geom_pointrange(aes(ymin=lower, ymax= upper), position=position_dodge(width=.9), size=1)+
  theme_bw() + 
  #facet_wrap(~dv, scales = "free_y") +
  ylab("Predicted Support") +
  xlab("Primary Measures") +
  geom_hline(aes(yintercept = yint), lty="dashed") +
  scale_shape_manual(name = "Included\nAncillary Measures", values = c(0:2, 5, 15, 17, 19, 3:4)) +
  scale_colour_manual(name = "Number of \nAncillary Measures", values = c("#5e3c99", "#e66101", "#b2abd2", "#fdb863")) +
  scale_linetype_manual(name = "Number of \nAncillary Measures", values = c("dotted", "dotdash", "dashed", "solid")) +
  theme(strip.background = element_rect(fill = 'white')) +
  theme(text = element_text(size=16)) +
  geom_vline(xintercept = 2.5, lty="dashed") +
  annotate("text", x = 2.4, y = .615, label = c('Command-and-Control'), hjust = 1)+
  annotate("text", x = 2.6, y = .615, label = c('Market-based'), hjust = 0)+
  guides(colour = guide_legend(override.aes = list(shape = NA)))

#Export Figure 3
ggsave(p_Anc_choice, filename = "./figs/packaging/JPP/JPP_Figure3.png", width = 12, height = 7, unit="in", dpi = 600)

##################### Appendix #################
####### Figures Appendix #######
#### Figure A3 ####
ggsave(p_sumAnc, filename = "./figs/packaging/JPP/JPP_FigureA3.png", width = 12, height = 7, unit="in", dpi = 600)

#### Figure A4 ####
p_Anc_rate <- ggplot(subset(fit, dv == "Rating"), aes(x= prim, y=fit, group = labels, shape = labels, colour = col, linetype = col))+
  geom_pointrange(aes(ymin=lower, ymax= upper), position=position_dodge(width=.9), size=1) +
  theme_bw() + 
  ylab("Predicted Support") +
  xlab("Primary Measures") +
  geom_hline(aes(yintercept = yint), lty="dashed") +
  scale_shape_manual(name = "Included\nAncillary Measures", values = c(0:2, 5, 15, 17, 19, 3:4)) +
  scale_colour_manual(name = "Number of \nAncillary Measures", values = c("#5e3c99", "#e66101", "#b2abd2", "#fdb863")) +
  scale_linetype_manual(name = "Number of \nAncillary Measures", values = c("dotted", "dotdash", "dashed", "solid")) +
  theme(strip.background = element_rect(fill = 'white')) +
  theme(text = element_text(size=16))+
  geom_vline(xintercept = 2.5, lty="dashed") +
  annotate("text", x = 2.4, y = 4.5, label = c('Command-and-Control'), hjust = 1)+
  annotate("text", x = 2.6, y = 4.5, label = c('Market-based'), hjust = 0) +
  guides(colour = guide_legend(override.aes = list(shape = NA)))

ggsave(p_Anc_rate, filename = "./figs/packaging/JPP/JPP_FigureA4.png", width = 12, height = 7, unit="in", dpi = 600)

#### Figure A5 and Figure A6 ####
# Margins ###

fit$relFit <- ifelse(fit$prim == "Registration\nBan" & fit$dv == "Choice", fit$fit - fit$fit[which(fit$prim == "Registration\nBan" & fit$labels == "Base Line" & fit$dv == "Choice")],
                     ifelse(fit$prim == "Registration\nBan" & fit$dv == "Rating", fit$fit - fit$fit[which(fit$prim == "Registration\nBan" & fit$labels == "Base Line" & fit$dv == "Rating")],
                            ifelse(fit$prim == "Emission\nDependent\nCar Tax" & fit$dv == "Choice", fit$fit - fit$fit[which(fit$prim == "Emission\nDependent\nCar Tax" & fit$labels == "Base Line" & fit$dv == "Choice")],
                                   ifelse(fit$prim == "Emission\nDependent\nCar Tax" & fit$dv == "Rating", fit$fit - fit$fit[which(fit$prim == "Emission\nDependent\nCar Tax" & fit$labels == "Base Line" & fit$dv == "Rating")],
                                          ifelse(fit$prim == "Stricter\nEnergy\nLabel" & fit$dv == "Choice", fit$fit - fit$fit[which(fit$prim == "Stricter\nEnergy\nLabel" & fit$labels == "Base Line" & fit$dv == "Choice")],
                                                 ifelse(fit$prim == "Stricter\nEnergy\nLabel" & fit$dv == "Rating", fit$fit - fit$fit[which(fit$prim == "Stricter\nEnergy\nLabel" & fit$labels == "Base Line" & fit$dv == "Rating")],
                                                        ifelse(fit$prim == "Road\nPricing" & fit$dv == "Choice", fit$fit - fit$fit[which(fit$prim == "Road\nPricing" & fit$labels == "Base Line" & fit$dv == "Choice")],
                                                               ifelse(fit$prim == "Road\nPricing" & fit$dv == "Rating", fit$fit - fit$fit[which(fit$prim == "Road\nPricing" & fit$labels == "Base Line" & fit$dv == "Rating")], NA))))))))

p_margins_Anc_rate <- ggplot(subset(fit, dv== "Rating"), aes(x= prim, y=relFit, group = labels, shape = labels, colour = col))+
  geom_pointrange(aes(ymin=relFit - 1.96*se, ymax= relFit + 1.96*se), position=position_dodge(width=.9), size=1)+
  theme_bw() + 
  ylab("Predicted Support") +
  xlab("Primary Measures") +
  geom_hline(yintercept = 0, lty="dashed") +
  scale_shape_manual(name = "Included\nAncillary Measures", values = c(0:2, 5, 15, 17, 19, 3:4)) +
  scale_colour_manual(name = "Number of \nAncillary Measures", values = c("#5e3c99", "#e66101", "#b2abd2", "#fdb863")) +
  theme(strip.background = element_rect(fill = 'white')) +
  theme(text = element_text(size=16)) +
  geom_vline(xintercept = 2.5, lty="dashed") +
  annotate("text", x = 2.4, y = .3, label = c('Command-and-Control'), hjust = 1)+
  annotate("text", x = 2.6, y = 0.3, label = c('Market-based'), hjust = 0)

p_margins_Anc_choice <- ggplot(subset(fit, dv== "Choice"), aes(x= prim, y=relFit, group = labels, shape = labels, colour = col))+
  geom_pointrange(aes(ymin=relFit - 1.96*se, ymax= relFit + 1.96*se), position=position_dodge(width=.9), size=1)+
  theme_bw() + 
  ylab("Predicted Support") +
  xlab("Primary Measures") +
  geom_hline(yintercept = 0, lty="dashed") +
  scale_shape_manual(name = "Included\nAncillary Measures", values = c(0:2, 5, 15, 17, 19, 3:4)) +
  scale_colour_manual(name = "Number of \nAncillary Measures", values = c("#5e3c99", "#e66101", "#b2abd2", "#fdb863")) +
  theme(strip.background = element_rect(fill = 'white')) +
  theme(text = element_text(size=16))+
  geom_vline(xintercept = 2.5, lty="dashed") +
  annotate("text", x = 2.4, y = .15, label = c('Command-and-Control'), hjust = 1)+
  annotate("text", x = 2.6, y = 0.15, label = c('Market-based'), hjust = 0)

ggsave(p_margins_Anc_choice, filename = "./figs/packaging/JPP/JPP_FigureA5.png", width = 12, height = 7, unit="in", dpi = 600)
ggsave(p_margins_Anc_rate, filename = "./figs/packaging/JPP/JPP_FigureA6.png", width = 12, height = 7, unit="in", dpi = 600)


#### Figure A1 ####
###################### number of primary and secondary measures in package

df_conj_out2$sumpm <- ifelse(df_conj_out2$Attrib1 == "_2", 1, 0) +
  ifelse(df_conj_out2$Attrib2 == "_2", 1, 0) +
  ifelse(df_conj_out2$Attrib3 == "_2", 1, 0) +
  ifelse(df_conj_out2$Attrib4 == "_2", 1, 0)

df_conj_out2$sumsm <- ifelse(df_conj_out2$Attrib6 == "_2", 1, 0) + 
  ifelse(df_conj_out2$Attrib8 == "_2", 1, 0) + 
  ifelse(df_conj_out2$Attrib9 == "_2", 1, 0)

df_conj_out2$sumsm_sq <- df_conj_out2$sumsm * df_conj_out2$sumsm
df_conj_out2$sumpm_sq <- df_conj_out2$sumpm * df_conj_out2$sumpm
df_conj_out2$sumpm_sumsm <- df_conj_out2$sumpm * df_conj_out2$sumsm

table(df_conj_out2$sumpm)

df_conj_out2$sum34 <-   ifelse(df_conj_out2$Attrib3 == "_2", 1, 0) +
  ifelse(df_conj_out2$Attrib4 == "_2", 1, 0)

table(df_conj_out2$sum34)

#rating visualisation
df_conj_out2$edu <- as.factor(df_conj_out2$edu)
m1_visual_rate <- lmer(rate ~ sumsm + sumpm + sumsm*sumpm + I(sumsm^2) + I(sumpm^2) + age+ gender+edu+usagecar + frame+ (1 | id), data = df_conj_out2)
summary(m1_visual_rate)

fit <- effect("sumsm", m1_visual_rate, xlevels = list(sumsm=c(0:3)), 
              x.var="sumsm", confidence.level = 0.95)

fit <- data.frame(fit)

ggplot(fit, aes(sumsm, fit)) +
  geom_line() +
  geom_line(aes(y=lower), lty="dashed") +
  geom_line(aes(y=upper), lty="dashed")

fit2 <- effect("sumpm", m1_visual_rate, xlevels = list(sumpm=c(0:4)), 
               x.var="sumpm", confidence.level = 0.95)

fit2 <- data.frame(fit2)

ggplot(fit2, aes(sumpm, fit)) +
  geom_line() +
  geom_line(aes(y=lower), lty="dashed") +
  geom_line(aes(y=upper), lty="dashed")

fit3 <- effect("sumsm:sumpm", m1_visual_rate, xlevels = list(sumpm=c(0:4),sumsm=c(0:3)), 
               x.var="sumpm", confidence.level = 0.95)

fit3 <- data.frame (fit3)

sumpm_sumsm <- ggplot(fit3, aes(x=sumpm, y=fit, colour=factor(sumsm))) +
  geom_pointrange(aes(min = lower, max = upper), position = position_dodge(width=.3) ) +
  theme_bw() + 
  scale_colour_manual(name = "Number of \nAncillary Measures", values = c("#5e3c99", "#e66101", "#b2abd2", "#fdb863")) +
  labs(x="Number of primary measures", y="Rating prediction", colour ="Number of ancillary measures")
sumpm_sumsm

ggsave("./figs/packaging/JPP/JPP_FigureA1.png", sumpm_sumsm, dpi = 600)

#### Figure A2 ####
sumsm_sumpm <- ggplot(fit3, aes(x=sumsm, y=fit, colour=factor(sumpm))) +
  geom_pointrange(aes(min = lower, max = upper), position = position_dodge(width=.3) ) +
  theme_bw() +
  scale_colour_manual(name = "Number of \nPrimary Measures", values = c("#000000", "#5e3c99", "#e66101", "#b2abd2", "#fdb863")) +
  labs(x="Number of ancillary measures", colour ="Number of primary measures")

ggsave("./figs/packaging/JPP/JPP_FigureA2.png", sumsm_sumpm, dpi = 600)

#### Interaction Primary Measures
#### Figure A7 ####

mprim_choice_4w <- lmer(choice ~ attrib1_lab * attrib2_lab * attrib3_lab * attrib4_lab + attrib5_lab + attrib6_lab + attrib7_lab + attrib8_lab + attrib9_lab + (1 | id), data = df_conj_out2)

fit <- effect("attrib1_lab:attrib2_lab:attrib3_lab:attrib4_lab", mprim_choice_4w, xlevels = list(attrib1_lab=levels(df_conj_out2$attrib1_lab),
                                                                                                 attrib2_lab=levels(df_conj_out2$attrib2_lab),
                                                                                                 attrib3_lab=levels(df_conj_out2$attrib3_lab),
                                                                                                 attrib4_lab=levels(df_conj_out2$attrib4_lab)), 
              x.var="attrib1_lab", confidence.level = 0.95)

fit <- data.frame(fit)

fit$nrattri <- factor(ifelse(fit$attrib1_lab == "2) Car Ban", 1, 0) + ifelse(fit$attrib3_lab == "2) Road Pricing", 1, 0) + ifelse(fit$attrib4_lab == "2) Car Tax", 1, 0) + ifelse(fit$attrib2_lab == "2) Tighter Energy Label", 1, 0))

fit$with <- factor(ifelse(fit$attrib2_lab == "1) No Tighter Energy Label" & fit$attrib1_lab == "1) No Car Ban" & fit$attrib4_lab == "1) No Car Tax", "Baseline",
                          ifelse(fit$attrib2_lab == "2) Tighter Energy Label" & fit$attrib1_lab == "1) No Car Ban" & fit$attrib4_lab == "1) No Car Tax", "w/ Tighter Energy Label",
                                 ifelse(fit$attrib2_lab == "1) No Tighter Energy Label" & fit$attrib1_lab == "2) Car Ban" & fit$attrib4_lab == "1) No Car Tax", "w/ Car Ban",
                                        ifelse(fit$attrib2_lab == "1) No Tighter Energy Label" & fit$attrib1_lab == "1) No Car Ban" & fit$attrib4_lab == "2) Car Tax", "w/ Car Tax",
                                               ifelse(fit$attrib2_lab == "2) Tighter Energy Label" & fit$attrib1_lab == "2) Car Ban" & fit$attrib4_lab == "1) No Car Tax", "w/ Car Ban & Tighter Energy Label",
                                                      ifelse(fit$attrib2_lab == "2) Tighter Energy Label" & fit$attrib1_lab == "1) No Car Ban" & fit$attrib4_lab == "2) Car Tax", "w/ Tighter Energy Label & Car Tax",
                                                             ifelse(fit$attrib2_lab == "1) No Tighter Energy Label" & fit$attrib1_lab == "2) Car Ban" & fit$attrib4_lab == "2) Car Tax", "w/ Car Ban & Car Tax",
                                                                    ifelse(fit$attrib2_lab == "2) Tighter Energy Label" & fit$attrib1_lab == "2) Car Ban" & fit$attrib4_lab == "2) Car Tax", "w/ Car Ban, Tighter Energy Label & Car Tax", NA)))))))),
                   levels = c("Baseline", "w/ Car Ban", "w/ Tighter Energy Label", "w/ Car Tax", "w/ Car Ban & Tighter Energy Label", "w/ Car Ban & Car Tax", "w/ Tighter Energy Label & Car Tax", "w/ Car Ban, Tighter Energy Label & Car Tax"),
                   labels = c("Baseline", "w/ Registration Ban", "w/ Stricter Energy Label", "w/ Emission Dependent\nCar Tax", "w/ Registration Ban & \nStricter Energy Label", "w/ Registration Ban &\nEmission Dependent\nCar Tax", "w/ Stricter Energy Label &\nEmission Dependent\nCar Tax", "w/ Registration Ban &\nStricter Energy Label &\nEmission Dependent\nCar Tax"))

p_prim_choice <- ggplot(subset(fit, nrattri != 0), aes(x= attrib3_lab, y=fit, group = with, colour = nrattri))+
  geom_pointrange(aes(ymin=lower, ymax= upper, shape = with), position=position_dodge(width=.5), size=1)+
  theme_bw() + 
  ylab("Predicted Support") +
  xlab("Road Pricing Included?") +
  geom_hline(yintercept = 0.5, lty="dashed") +
  scale_colour_manual(name = "Number of\n Primary Measures", values = c("#5e3c99", "#e66101", "#b2abd2", "#fdb863")) +
  scale_shape_manual(name = "Included\nPrimary Measures", values = c(0:2, 5, 15, 17, 19, 3:4)) +
  theme(strip.background = element_rect(fill = 'white')) +
  theme(text = element_text(size=16)) +
  scale_x_discrete(labels = c("No Road Pricing", "Road Pricing")) +
  NULL
p_prim_choice

ggsave("./figs/packaging/JPP/JPP_FigureA7.png", p_prim_choice, dpi = 600)

#### Figure A8 ####
mprim_rate_4w <- lmer(rate ~ attrib1_lab * attrib2_lab * attrib3_lab * attrib4_lab + attrib5_lab + attrib6_lab + attrib7_lab + attrib8_lab + attrib9_lab + (1 | id), data = df_conj_out2)

fit <- effect("attrib1_lab:attrib2_lab:attrib3_lab:attrib4_lab", mprim_rate_4w, xlevels = list(attrib1_lab=levels(df_conj_out2$attrib1_lab),
                                                                                               attrib2_lab=levels(df_conj_out2$attrib2_lab),
                                                                                               attrib3_lab=levels(df_conj_out2$attrib3_lab),
                                                                                               attrib4_lab=levels(df_conj_out2$attrib4_lab)), 
              x.var="attrib1_lab", confidence.level = 0.95)

fit <- data.frame(fit)

fit$nrattri <- factor(ifelse(fit$attrib1_lab == "2) Car Ban", 1, 0) + ifelse(fit$attrib3_lab == "2) Road Pricing", 1, 0) + ifelse(fit$attrib4_lab == "2) Car Tax", 1, 0) + ifelse(fit$attrib2_lab == "2) Tighter Energy Label", 1, 0))

fit$with <- factor(ifelse(fit$attrib2_lab == "1) No Tighter Energy Label" & fit$attrib1_lab == "1) No Car Ban" & fit$attrib4_lab == "1) No Car Tax", "Baseline",
                          ifelse(fit$attrib2_lab == "2) Tighter Energy Label" & fit$attrib1_lab == "1) No Car Ban" & fit$attrib4_lab == "1) No Car Tax", "w/ Tighter Energy Label",
                                 ifelse(fit$attrib2_lab == "1) No Tighter Energy Label" & fit$attrib1_lab == "2) Car Ban" & fit$attrib4_lab == "1) No Car Tax", "w/ Car Ban",
                                        ifelse(fit$attrib2_lab == "1) No Tighter Energy Label" & fit$attrib1_lab == "1) No Car Ban" & fit$attrib4_lab == "2) Car Tax", "w/ Car Tax",
                                               ifelse(fit$attrib2_lab == "2) Tighter Energy Label" & fit$attrib1_lab == "2) Car Ban" & fit$attrib4_lab == "1) No Car Tax", "w/ Car Ban & Tighter Energy Label",
                                                      ifelse(fit$attrib2_lab == "2) Tighter Energy Label" & fit$attrib1_lab == "1) No Car Ban" & fit$attrib4_lab == "2) Car Tax", "w/ Tighter Energy Label & Car Tax",
                                                             ifelse(fit$attrib2_lab == "1) No Tighter Energy Label" & fit$attrib1_lab == "2) Car Ban" & fit$attrib4_lab == "2) Car Tax", "w/ Car Ban & Car Tax",
                                                                    ifelse(fit$attrib2_lab == "2) Tighter Energy Label" & fit$attrib1_lab == "2) Car Ban" & fit$attrib4_lab == "2) Car Tax", "w/ Car Ban, Tighter Energy Label & Car Tax", NA)))))))),
                   levels = c("Baseline", "w/ Car Ban", "w/ Tighter Energy Label", "w/ Car Tax", "w/ Car Ban & Tighter Energy Label", "w/ Car Ban & Car Tax", "w/ Tighter Energy Label & Car Tax", "w/ Car Ban, Tighter Energy Label & Car Tax"),
                   labels = c("Baseline", "w/ Registration Ban", "w/ Stricter Energy Label", "w/ Emission Dependent\nCar Tax", "w/ Registration Ban & \nStricter Energy Label", "w/ Registration Ban &\nEmission Dependent\nCar Tax", "w/ Stricter Energy Label &\nEmission Dependent\nCar Tax", "w/ Registration Ban &\nStricter Energy Label &\nEmission Dependent\nCar Tax"))




p_prim_rate <- ggplot(subset(fit, nrattri != 0), aes(x= attrib3_lab, y=fit, group = with, colour = nrattri))+
  geom_pointrange(aes(ymin=lower, ymax= upper, shape = with), position=position_dodge(width=.5), size=1)+
  theme_bw() + 
  ylab("Predicted Support") +
  xlab("Road Pricing Included?") +
  geom_hline(yintercept = 4, lty="dashed") +
  scale_colour_manual(name = "Number of\n Primary Measures", values = c("#5e3c99", "#e66101", "#b2abd2", "#fdb863")) +
  scale_shape_manual(name = "Included\nPrimary Measures", values = c(0:2, 5, 15, 17, 19, 3:4)) +
  theme(strip.background = element_rect(fill = 'white')) +
  theme(text = element_text(size=16)) +
  scale_x_discrete(labels = c("No Road Pricing", "Road Pricing")) +
  NULL
p_prim_rate

ggsave("./figs/packaging/JPP/JPP_FigureA8.png", p_prim_rate, dpi = 600)

######## Figure A9 and A10 ####### 
#Learning Effects by round #
######## Choice
#round 1
df_conj_out2r1 <- subset(df_conj_out2, round == 1)

c_base1 <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r1, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(c_base1)

c_attribute1 <- c(0, c_base1$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, c_base1$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, c_base1$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, c_base1$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, c_base1$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, c_base1$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, c_base1$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, c_base1$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, c_base1$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base1$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], c_base1$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, c_base1$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base1$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], c_base1$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, c_base1$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, c_base1$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, c_base1$estimates$attrib7lab["AMCE","attrib7lab22030"], c_base1$estimates$attrib7lab["AMCE","attrib7lab32035"], c_base1$estimates$attrib7lab["AMCE","attrib7lab42040"], c_base1$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, c_base1$estimates$attrib7lab["Std. Error","attrib7lab22030"], c_base1$estimates$attrib7lab["Std. Error","attrib7lab32035"], c_base1$estimates$attrib7lab["Std. Error","attrib7lab42040"], c_base1$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, c_base1$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, c_base1$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, c_base1$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, c_base1$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_c_base1 <- data.frame(c,s)
df_c_base1$att <- c(c_base1$attributes$attrib1lab[1], c_base1$attributes$attrib1lab[2], c_base1$attributes$attrib2lab[1], c_base1$attributes$attrib2lab[2], c_base1$attributes$attrib3lab[1], c_base1$attributes$attrib3lab[2], c_base1$attributes$attrib4lab[1], c_base1$attributes$attrib4lab[2], c_base1$attributes$attrib5lab[1], c_base1$attributes$attrib5lab[2:4], c_base1$attributes$attrib6lab[1], c_base1$attributes$attrib6lab[2], c_base1$attributes$attrib7lab[1], c_base1$attributes$attrib7lab[2:5], c_base1$attributes$attrib8lab[1], c_base1$attributes$attrib8lab[2], c_base1$attributes$attrib9lab[1], c_base1$attributes$attrib9lab[2])
df_c_base1$att <- substr(df_c_base1$att, 2, stop = nchar(df_c_base1$att))
df_c_base1$att <- factor(df_c_base1$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))

#round 2
df_conj_out2r2 <- subset(df_conj_out2, round == 2)

c_base2 <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r2, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(c_base2)

c_attribute1 <- c(0, c_base2$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, c_base2$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, c_base2$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, c_base2$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, c_base2$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, c_base2$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, c_base2$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, c_base2$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, c_base2$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base2$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], c_base2$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, c_base2$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base2$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], c_base2$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, c_base2$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, c_base2$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, c_base2$estimates$attrib7lab["AMCE","attrib7lab22030"], c_base2$estimates$attrib7lab["AMCE","attrib7lab32035"], c_base2$estimates$attrib7lab["AMCE","attrib7lab42040"], c_base2$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, c_base2$estimates$attrib7lab["Std. Error","attrib7lab22030"], c_base2$estimates$attrib7lab["Std. Error","attrib7lab32035"], c_base2$estimates$attrib7lab["Std. Error","attrib7lab42040"], c_base2$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, c_base2$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, c_base2$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, c_base2$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, c_base2$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_c_base2 <- data.frame(c,s)
df_c_base2$att <- c(c_base2$attributes$attrib1lab[1], c_base2$attributes$attrib1lab[2], c_base2$attributes$attrib2lab[1], c_base2$attributes$attrib2lab[2], c_base2$attributes$attrib3lab[1], c_base2$attributes$attrib3lab[2], c_base2$attributes$attrib4lab[1], c_base2$attributes$attrib4lab[2], c_base2$attributes$attrib5lab[1], c_base2$attributes$attrib5lab[2:4], c_base2$attributes$attrib6lab[1], c_base2$attributes$attrib6lab[2], c_base2$attributes$attrib7lab[1], c_base2$attributes$attrib7lab[2:5], c_base2$attributes$attrib8lab[1], c_base2$attributes$attrib8lab[2], c_base2$attributes$attrib9lab[1], c_base2$attributes$attrib9lab[2])
df_c_base2$att <- substr(df_c_base2$att, 2, stop = nchar(df_c_base2$att))
df_c_base2$att <- factor(df_c_base2$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))

#round 3
df_conj_out2r3 <- subset(df_conj_out2, round == 3)

c_base3 <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r3, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(c_base3)

c_attribute1 <- c(0, c_base3$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, c_base3$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, c_base3$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, c_base3$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, c_base3$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, c_base3$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, c_base3$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, c_base3$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, c_base3$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base3$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], c_base3$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, c_base3$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base3$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], c_base3$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, c_base3$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, c_base3$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, c_base3$estimates$attrib7lab["AMCE","attrib7lab22030"], c_base3$estimates$attrib7lab["AMCE","attrib7lab32035"], c_base3$estimates$attrib7lab["AMCE","attrib7lab42040"], c_base3$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, c_base3$estimates$attrib7lab["Std. Error","attrib7lab22030"], c_base3$estimates$attrib7lab["Std. Error","attrib7lab32035"], c_base3$estimates$attrib7lab["Std. Error","attrib7lab42040"], c_base3$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, c_base3$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, c_base3$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, c_base3$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, c_base3$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_c_base3 <- data.frame(c,s)
df_c_base3$att <- c(c_base3$attributes$attrib1lab[1], c_base3$attributes$attrib1lab[2], c_base3$attributes$attrib2lab[1], c_base3$attributes$attrib2lab[2], c_base3$attributes$attrib3lab[1], c_base3$attributes$attrib3lab[2], c_base3$attributes$attrib4lab[1], c_base3$attributes$attrib4lab[2], c_base3$attributes$attrib5lab[1], c_base3$attributes$attrib5lab[2:4], c_base3$attributes$attrib6lab[1], c_base3$attributes$attrib6lab[2], c_base3$attributes$attrib7lab[1], c_base3$attributes$attrib7lab[2:5], c_base3$attributes$attrib8lab[1], c_base3$attributes$attrib8lab[2], c_base3$attributes$attrib9lab[1], c_base3$attributes$attrib9lab[2])
df_c_base3$att <- substr(df_c_base3$att, 2, stop = nchar(df_c_base3$att))
df_c_base3$att <- factor(df_c_base3$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))

#round 4
df_conj_out2r4 <- subset(df_conj_out2, round == 4)

c_base4 <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r4, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(c_base4)

c_attribute1 <- c(0, c_base4$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, c_base4$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, c_base4$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, c_base4$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, c_base4$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, c_base4$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, c_base4$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, c_base4$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, c_base4$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base4$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], c_base4$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, c_base4$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base4$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], c_base4$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, c_base4$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, c_base4$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, c_base4$estimates$attrib7lab["AMCE","attrib7lab22030"], c_base4$estimates$attrib7lab["AMCE","attrib7lab32035"], c_base4$estimates$attrib7lab["AMCE","attrib7lab42040"], c_base4$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, c_base4$estimates$attrib7lab["Std. Error","attrib7lab22030"], c_base4$estimates$attrib7lab["Std. Error","attrib7lab32035"], c_base4$estimates$attrib7lab["Std. Error","attrib7lab42040"], c_base4$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, c_base4$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, c_base4$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, c_base4$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, c_base4$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_c_base4 <- data.frame(c,s)
df_c_base4$att <- c(c_base4$attributes$attrib1lab[1], c_base4$attributes$attrib1lab[2], c_base4$attributes$attrib2lab[1], c_base4$attributes$attrib2lab[2], c_base4$attributes$attrib3lab[1], c_base4$attributes$attrib3lab[2], c_base4$attributes$attrib4lab[1], c_base4$attributes$attrib4lab[2], c_base4$attributes$attrib5lab[1], c_base4$attributes$attrib5lab[2:4], c_base4$attributes$attrib6lab[1], c_base4$attributes$attrib6lab[2], c_base4$attributes$attrib7lab[1], c_base4$attributes$attrib7lab[2:5], c_base4$attributes$attrib8lab[1], c_base4$attributes$attrib8lab[2], c_base4$attributes$attrib9lab[1], c_base4$attributes$attrib9lab[2])
df_c_base4$att <- substr(df_c_base4$att, 2, stop = nchar(df_c_base4$att))
df_c_base4$att <- factor(df_c_base4$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))

#round 5
df_conj_out2r5 <- subset(df_conj_out2, round == 5)

c_base5 <- amce(choice ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r5, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(c_base5)

c_attribute1 <- c(0, c_base5$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, c_base5$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, c_base5$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, c_base5$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, c_base5$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, c_base5$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, c_base5$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, c_base5$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, c_base5$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base5$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], c_base5$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, c_base5$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], c_base5$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], c_base5$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, c_base5$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, c_base5$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, c_base5$estimates$attrib7lab["AMCE","attrib7lab22030"], c_base5$estimates$attrib7lab["AMCE","attrib7lab32035"], c_base5$estimates$attrib7lab["AMCE","attrib7lab42040"], c_base5$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, c_base5$estimates$attrib7lab["Std. Error","attrib7lab22030"], c_base5$estimates$attrib7lab["Std. Error","attrib7lab32035"], c_base5$estimates$attrib7lab["Std. Error","attrib7lab42040"], c_base5$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, c_base5$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, c_base5$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, c_base5$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, c_base5$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_c_base5 <- data.frame(c,s)
df_c_base5$att <- c(c_base5$attributes$attrib1lab[1], c_base5$attributes$attrib1lab[2], c_base5$attributes$attrib2lab[1], c_base5$attributes$attrib2lab[2], c_base5$attributes$attrib3lab[1], c_base5$attributes$attrib3lab[2], c_base5$attributes$attrib4lab[1], c_base5$attributes$attrib4lab[2], c_base5$attributes$attrib5lab[1], c_base5$attributes$attrib5lab[2:4], c_base5$attributes$attrib6lab[1], c_base5$attributes$attrib6lab[2], c_base5$attributes$attrib7lab[1], c_base5$attributes$attrib7lab[2:5], c_base5$attributes$attrib8lab[1], c_base5$attributes$attrib8lab[2], c_base5$attributes$attrib9lab[1], c_base5$attributes$attrib9lab[2])
df_c_base5$att <- substr(df_c_base5$att, 2, stop = nchar(df_c_base5$att))
df_c_base5$att <- factor(df_c_base5$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))


df_rounds <- rbind(df_c_base1, df_c_base2, df_c_base3, df_c_base4, df_c_base5)
df_rounds$rnd <- rep(c("Round 1", "Round 2", "Round 3", "Round 4", "Round 5"), each=nrow(df_c_base1))

rounds_output <- ggplot(df_rounds, aes(x=att, y = c, group = rnd, shape = rnd)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(legend.position = "bottom") +
  coord_flip() +
  ylab("Change in E[Y]") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") + 
  #facet_wrap(~rnd, scales = "free_x") +
  scale_shape_manual(name="Round", values = c(0:4)) + 
  theme(strip.background = element_rect(fill = 'white'))

########Figure A9
ggsave("./figs/packaging/JPP/JPP_FigureA9.png", rounds_output, width = 9, height = 11, dpi = 600)

###### Rating
r_base1 <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r1, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(r_base1)

c_attribute1 <- c(0, r_base1$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, r_base1$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, r_base1$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, r_base1$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, r_base1$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, r_base1$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, r_base1$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, r_base1$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, r_base1$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base1$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], r_base1$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, r_base1$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base1$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], r_base1$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, r_base1$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, r_base1$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, r_base1$estimates$attrib7lab["AMCE","attrib7lab22030"], r_base1$estimates$attrib7lab["AMCE","attrib7lab32035"], r_base1$estimates$attrib7lab["AMCE","attrib7lab42040"], r_base1$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, r_base1$estimates$attrib7lab["Std. Error","attrib7lab22030"], r_base1$estimates$attrib7lab["Std. Error","attrib7lab32035"], r_base1$estimates$attrib7lab["Std. Error","attrib7lab42040"], r_base1$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, r_base1$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, r_base1$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, r_base1$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, r_base1$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_r_base1 <- data.frame(c,s)
df_r_base1$att <- c(r_base1$attributes$attrib1lab[1], r_base1$attributes$attrib1lab[2], r_base1$attributes$attrib2lab[1], r_base1$attributes$attrib2lab[2], r_base1$attributes$attrib3lab[1], r_base1$attributes$attrib3lab[2], r_base1$attributes$attrib4lab[1], r_base1$attributes$attrib4lab[2], r_base1$attributes$attrib5lab[1], r_base1$attributes$attrib5lab[2:4], r_base1$attributes$attrib6lab[1], r_base1$attributes$attrib6lab[2], r_base1$attributes$attrib7lab[1], r_base1$attributes$attrib7lab[2:5], r_base1$attributes$attrib8lab[1], r_base1$attributes$attrib8lab[2], r_base1$attributes$attrib9lab[1], r_base1$attributes$attrib9lab[2])
df_r_base1$att <- substr(df_r_base1$att, 2, stop = nchar(df_r_base1$att))
df_r_base1$att <- factor(df_r_base1$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))

#round 2
r_base2 <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r2, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(r_base2)

c_attribute1 <- c(0, r_base2$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, r_base2$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, r_base2$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, r_base2$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, r_base2$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, r_base2$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, r_base2$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, r_base2$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, r_base2$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base2$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], r_base2$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, r_base2$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base2$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], r_base2$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, r_base2$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, r_base2$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, r_base2$estimates$attrib7lab["AMCE","attrib7lab22030"], r_base2$estimates$attrib7lab["AMCE","attrib7lab32035"], r_base2$estimates$attrib7lab["AMCE","attrib7lab42040"], r_base2$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, r_base2$estimates$attrib7lab["Std. Error","attrib7lab22030"], r_base2$estimates$attrib7lab["Std. Error","attrib7lab32035"], r_base2$estimates$attrib7lab["Std. Error","attrib7lab42040"], r_base2$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, r_base2$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, r_base2$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, r_base2$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, r_base2$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_r_base2 <- data.frame(c,s)
df_r_base2$att <- c(r_base2$attributes$attrib1lab[1], r_base2$attributes$attrib1lab[2], r_base2$attributes$attrib2lab[1], r_base2$attributes$attrib2lab[2], r_base2$attributes$attrib3lab[1], r_base2$attributes$attrib3lab[2], r_base2$attributes$attrib4lab[1], r_base2$attributes$attrib4lab[2], r_base2$attributes$attrib5lab[1], r_base2$attributes$attrib5lab[2:4], r_base2$attributes$attrib6lab[1], r_base2$attributes$attrib6lab[2], r_base2$attributes$attrib7lab[1], r_base2$attributes$attrib7lab[2:5], r_base2$attributes$attrib8lab[1], r_base2$attributes$attrib8lab[2], r_base2$attributes$attrib9lab[1], r_base2$attributes$attrib9lab[2])
df_r_base2$att <- substr(df_r_base2$att, 2, stop = nchar(df_r_base2$att))
df_r_base2$att <- factor(df_r_base2$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))  
#round 3
r_base3 <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r3, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(r_base3)

c_attribute1 <- c(0, r_base3$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, r_base3$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, r_base3$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, r_base3$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, r_base3$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, r_base3$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, r_base3$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, r_base3$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, r_base3$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base3$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], r_base3$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, r_base3$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base3$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], r_base3$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, r_base3$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, r_base3$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, r_base3$estimates$attrib7lab["AMCE","attrib7lab22030"], r_base3$estimates$attrib7lab["AMCE","attrib7lab32035"], r_base3$estimates$attrib7lab["AMCE","attrib7lab42040"], r_base3$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, r_base3$estimates$attrib7lab["Std. Error","attrib7lab22030"], r_base3$estimates$attrib7lab["Std. Error","attrib7lab32035"], r_base3$estimates$attrib7lab["Std. Error","attrib7lab42040"], r_base3$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, r_base3$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, r_base3$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, r_base3$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, r_base3$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_r_base3 <- data.frame(c,s)
df_r_base3$att <- c(r_base3$attributes$attrib1lab[1], r_base3$attributes$attrib1lab[2], r_base3$attributes$attrib2lab[1], r_base3$attributes$attrib2lab[2], r_base3$attributes$attrib3lab[1], r_base3$attributes$attrib3lab[2], r_base3$attributes$attrib4lab[1], r_base3$attributes$attrib4lab[2], r_base3$attributes$attrib5lab[1], r_base3$attributes$attrib5lab[2:4], r_base3$attributes$attrib6lab[1], r_base3$attributes$attrib6lab[2], r_base3$attributes$attrib7lab[1], r_base3$attributes$attrib7lab[2:5], r_base3$attributes$attrib8lab[1], r_base3$attributes$attrib8lab[2], r_base3$attributes$attrib9lab[1], r_base3$attributes$attrib9lab[2])
df_r_base3$att <- substr(df_r_base3$att, 2, stop = nchar(df_r_base3$att))
df_r_base3$att <- factor(df_r_base3$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))  
#round 4
r_base4 <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r4, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(r_base4)

c_attribute1 <- c(0, r_base4$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, r_base4$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, r_base4$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, r_base4$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, r_base4$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, r_base4$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, r_base4$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, r_base4$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, r_base4$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base4$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], r_base4$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, r_base4$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base4$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], r_base4$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, r_base4$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, r_base4$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, r_base4$estimates$attrib7lab["AMCE","attrib7lab22030"], r_base4$estimates$attrib7lab["AMCE","attrib7lab32035"], r_base4$estimates$attrib7lab["AMCE","attrib7lab42040"], r_base4$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, r_base4$estimates$attrib7lab["Std. Error","attrib7lab22030"], r_base4$estimates$attrib7lab["Std. Error","attrib7lab32035"], r_base4$estimates$attrib7lab["Std. Error","attrib7lab42040"], r_base4$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, r_base4$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, r_base4$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, r_base4$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, r_base4$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_r_base4 <- data.frame(c,s)
df_r_base4$att <- c(r_base4$attributes$attrib1lab[1], r_base4$attributes$attrib1lab[2], r_base4$attributes$attrib2lab[1], r_base4$attributes$attrib2lab[2], r_base4$attributes$attrib3lab[1], r_base4$attributes$attrib3lab[2], r_base4$attributes$attrib4lab[1], r_base4$attributes$attrib4lab[2], r_base4$attributes$attrib5lab[1], r_base4$attributes$attrib5lab[2:4], r_base4$attributes$attrib6lab[1], r_base4$attributes$attrib6lab[2], r_base4$attributes$attrib7lab[1], r_base4$attributes$attrib7lab[2:5], r_base4$attributes$attrib8lab[1], r_base4$attributes$attrib8lab[2], r_base4$attributes$attrib9lab[1], r_base4$attributes$attrib9lab[2])
df_r_base4$att <- substr(df_r_base4$att, 2, stop = nchar(df_r_base4$att))
df_r_base4$att <- factor(df_r_base4$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))  
#round 5
r_base5 <- amce(rate ~ attrib1_lab +  attrib2_lab +  attrib3_lab +  attrib4_lab +  
                  attrib5_lab +  attrib6_lab +  attrib7_lab +  attrib8_lab +  
                  attrib9_lab, data = df_conj_out2r5, 
                cluster=TRUE, respondent.id="id", 
                design="uniform")

plot(r_base5)

c_attribute1 <- c(0, r_base5$estimates$attrib1lab["AMCE","attrib1lab2CarBan"])
s_attribute1 <- c(0, r_base5$estimates$attrib1lab["Std. Error","attrib1lab2CarBan"])
c_attribute2 <- c(0, r_base5$estimates$attrib2lab["AMCE","attrib2lab2TighterEnergyLabel"])
s_attribute2 <- c(0, r_base5$estimates$attrib2lab["Std. Error","attrib2lab2TighterEnergyLabel"])
c_attribute3 <- c(0, r_base5$estimates$attrib3lab["AMCE","attrib3lab2RoadPricing"])
s_attribute3 <- c(0, r_base5$estimates$attrib3lab["Std. Error","attrib3lab2RoadPricing"])
c_attribute4 <- c(0, r_base5$estimates$attrib4lab["AMCE","attrib4lab2CarTax"])
s_attribute4 <- c(0, r_base5$estimates$attrib4lab["Std. Error","attrib4lab2CarTax"])
c_attribute5 <- c(0, r_base5$estimates$attrib5lab["AMCE","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base5$estimates$attrib5lab["AMCE","attrib5lab3GeneralBudget"], r_base5$estimates$attrib5lab["AMCE","attrib5labPublicTransportationPriceReduction"])
s_attribute5 <- c(0, r_base5$estimates$attrib5lab["Std. Error","attrib5lab2StreetBuildingMaintenanceandPublicTransportInfrastructure"], r_base5$estimates$attrib5lab["Std. Error","attrib5lab3GeneralBudget"], r_base5$estimates$attrib5lab["Std. Error","attrib5labPublicTransportationPriceReduction"])
c_attribute6 <- c(0, r_base5$estimates$attrib6lab["AMCE","attrib6lab2Campaign"])
s_attribute6 <- c(0, r_base5$estimates$attrib6lab["Std. Error","attrib6lab2Campaign"])
c_attribute7 <- c(0, r_base5$estimates$attrib7lab["AMCE","attrib7lab22030"], r_base5$estimates$attrib7lab["AMCE","attrib7lab32035"], r_base5$estimates$attrib7lab["AMCE","attrib7lab42040"], r_base5$estimates$attrib7lab["AMCE","attrib7lab52045"])
s_attribute7 <- c(0, r_base5$estimates$attrib7lab["Std. Error","attrib7lab22030"], r_base5$estimates$attrib7lab["Std. Error","attrib7lab32035"], r_base5$estimates$attrib7lab["Std. Error","attrib7lab42040"], r_base5$estimates$attrib7lab["Std. Error","attrib7lab52045"])
c_attribute8 <- c(0, r_base5$estimates$attrib8lab["AMCE","attrib8lab2TestPhase"])
s_attribute8 <- c(0, r_base5$estimates$attrib8lab["Std. Error","attrib8lab2TestPhase"])
c_attribute9 <- c(0, r_base5$estimates$attrib9lab["AMCE","attrib9lab2FederalGovtandExperts"])
s_attribute9 <- c(0, r_base5$estimates$attrib9lab["Std. Error","attrib9lab2FederalGovtandExperts"])

c <- c(c_attribute1, c_attribute2, c_attribute3, c_attribute4, c_attribute5, c_attribute6, c_attribute7, c_attribute8, c_attribute9)
s <- c(s_attribute1, s_attribute2, s_attribute3, s_attribute4, s_attribute5, s_attribute6, s_attribute7, s_attribute8, s_attribute9)

df_r_base5 <- data.frame(c,s)
df_r_base5$att <- c(r_base5$attributes$attrib1lab[1], r_base5$attributes$attrib1lab[2], r_base5$attributes$attrib2lab[1], r_base5$attributes$attrib2lab[2], r_base5$attributes$attrib3lab[1], r_base5$attributes$attrib3lab[2], r_base5$attributes$attrib4lab[1], r_base5$attributes$attrib4lab[2], r_base5$attributes$attrib5lab[1], r_base5$attributes$attrib5lab[2:4], r_base5$attributes$attrib6lab[1], r_base5$attributes$attrib6lab[2], r_base5$attributes$attrib7lab[1], r_base5$attributes$attrib7lab[2:5], r_base5$attributes$attrib8lab[1], r_base5$attributes$attrib8lab[2], r_base5$attributes$attrib9lab[1], r_base5$attributes$attrib9lab[2])
df_r_base5$att <- substr(df_r_base5$att, 2, stop = nchar(df_r_base5$att))
df_r_base5$att <- factor(df_r_base5$att, levels = c(rev(c("NoCarBan", "CarBan", "NoTighterEnergyLabel", "TighterEnergyLabel", "NoRoadPricing", "RoadPricing", "NoCarTax", "CarTax", "StreetBuildingandMaintenance", "StreetBuildingMaintenanceandPublicTransportInfrastructure", "GeneralBudget", "ublicTransportationPriceReduction", "NoCampaign", "Campaign", "2025", "2030", "2035", "2040", "2045", "NoTestPhase", "TestPhase", "FedealGovt", "FederalGovtandExperts"))),
                         labels = c(rev(c("Baseline: No Registration Ban", "Registration Ban", "Baseline: No Stricter Energy Label", "Stricter Energy Label", "Baseline: No Road Pricing", "Road Pricing", "Baseline: No Emission Dependent Car Tax", "Emission Dependent Car Tax", "Baseline: Street Building and Maintenance", "Road and PT Infrastructure", "General Budget", "PT Price Reduction", "Baseline: No Campaign", "Campaign", "Baseline: 2025", "2030", "2035", "2040", "2045", "Baseline: No Test Phase", "Test Phase", "Baseline: No Expert Involvement", "Expert Involvement"))))  

df_rounds_r <- rbind(df_r_base1, df_r_base2, df_r_base3, df_r_base4, df_r_base5)
df_rounds_r$rnd <- rep(c("Round 1", "Round 2", "Round 3", "Round 4", "Round 5"), each=nrow(df_r_base1))

rounds_output_r <- ggplot(df_rounds_r, aes(x=att, y = c, group = rnd, shape = rnd)) +
  geom_pointrange(aes(min = c - 1.95 * s, max = c + 1.95 * s), position = position_dodge(width=.5)) +
  theme_bw() + 
  theme(legend.position = "bottom") +
  coord_flip() +
  ylab("Change in E[Y]") +
  xlab("Attributes") +
  geom_hline(yintercept = 0, lty="dashed") + 
  scale_shape_manual(name="Round", values = c(0:4)) + 
  theme(strip.background = element_rect(fill = 'white'))

#####FigureA10
ggsave("./figs/packaging/JPP/JPP_FigureA10.png", rounds_output_r, width = 9, height = 11, dpi = 600)


####### Tables Appendix ########
############ Table A1 #############
mCB_rateI2 <- lmer(rate ~ attrib1_lab *attrib6_lab + attrib1_lab* attrib8_lab + attrib1_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mCB_rateI3 <- lmer(rate ~ attrib1_lab *attrib6_lab * attrib8_lab +
                     attrib1_lab * attrib6_lab * attrib9_lab +
                     attrib1_lab * attrib8_lab * attrib9_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mCB_rateI4 <- lmer(rate ~ attrib1_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

htmlreg(list(mCB_rateI2, mCB_rateI3, mCB_rateI4),
        file = "./tables/JPP/A1_PackagingCB_Rating.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Car Ban", "Campaign", "Expert Involvement",
                              "Test Phase", "Energy Label", "Road Pricing", "Car Tax",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Car Ban x Campaign", "Car Ban x Expert Involvement",
                              "Car Ban x Test Phase", "Campaign x Test Phase", "Campaign x Expert Involvement",
                              "Test Phase x Expert Involvement", "Car Ban x Campaign x Test Phase",
                              "Car Ban x Campaign x Expert Involvement", "Car Ban x Test Phase x Expert Involvement",
                              "Campaign x Test Phase x Expert Involvement", "Car Ban x Campaign x Test Phase x Expert Involvement"),
        custom.model.names = c("Two-Way Interaction", "Three Way Interaction", "Four Way Interaction"),
        stars = c(0.01, 0.05, 0.1))

############ Table A2 #############
mCT_rateI2 <- lmer(rate ~ attrib4_lab *attrib6_lab + attrib4_lab* attrib8_lab + attrib4_lab* attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mCT_rateI3 <- lmer(rate ~ attrib4_lab *attrib6_lab * attrib8_lab+
                     attrib4_lab *attrib6_lab * attrib9_lab +
                     attrib4_lab *attrib9_lab * attrib8_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mCT_rateI4 <- lmer(rate ~ attrib4_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

htmlreg(list(mCT_rateI2, mCT_rateI3, mCT_rateI4),
        file = "./tables/JPP/A2_PackagingCT_Rating.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Car Tax", "Campaign", "Expert Involvement",
                              "Test Phase", "Car Ban", "Energy Label", "Road Pricing",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Car Tax x Campaign", "Car Tax x Expert Involvement",
                              "Car Tax x Test Phase", "Campaign x Test Phase", "Campaign x Expert Involvement",
                              "Test Phase x Expert Involvement", "Car Tax x Campaign x Test Phase",
                              "Car Tax x Campaign x Expert Involvement", "Car Tax x Test Phase x Expert Involvement",
                              "Campaign x Test Phase x Expert Involvement", "Car Tax x Campaign x Test Phase x Expert Involvement"),
        custom.model.names = c("Two-Way Interaction", "Three Way Interaction", "Four Way Interaction"),
        stars = c(0.01, 0.05, 0.1))

############ Table A3 #############
mEL_rateI2 <- lmer(rate ~ attrib2_lab *attrib6_lab + attrib2_lab * attrib8_lab + attrib2_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mEL_rateI3 <- lmer(rate ~ attrib2_lab *attrib6_lab * attrib8_lab +
                     attrib2_lab *attrib6_lab * attrib9_lab +
                     attrib2_lab *attrib9_lab * attrib8_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mEL_rateI4 <- lmer(rate ~ attrib2_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

htmlreg(list(mEL_rateI2, mEL_rateI3, mEL_rateI4),
        file = "./tables/JPP/A3_PackagingEL_Rating.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Energy Label", "Campaign", "Expert Involvement",
                              "Test Phase", "Car Ban", "Road Pricing", "Car Tax",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Energy Label x Campaign", "Energy Label x Expert Involvement",
                              "Energy Label x Test Phase", "Campaign x Test Phase", "Campaign x Expert Involvement",
                              "Test Phase x Expert Involvement", "Energy Label x Campaign x Test Phase",
                              "Energy Label x Campaign x Expert Involvement", "Energy Label x Test Phase x Expert Involvement",
                              "Campaign x Test Phase x Expert Involvement", "Energy Label x Campaign x Test Phase x Expert Involvement"),
        custom.model.names = c("Two-Way Interaction", "Three Way Interaction", "Four Way Interaction"),
        stars = c(0.01, 0.05, 0.1))

############ Table A4 #############
mRP_rateI2 <- lmer(rate ~ attrib3_lab *attrib6_lab + attrib3_lab * attrib9_lab + attrib3_lab* attrib8_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mRP_rateI3 <- lmer(rate ~ attrib3_lab *attrib6_lab * attrib8_lab+ attrib3_lab *attrib6_lab * attrib9_lab + attrib3_lab* attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mRP_rateI4 <- lmer(rate ~ attrib3_lab *attrib6_lab * attrib8_lab*  attrib9_lab  + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

htmlreg(list(mRP_rateI2, mRP_rateI3, mRP_rateI4),
        file = "./tables/JPP/A4_PackagingRP_Rating.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Road Pricing", "Campaign", "Expert Involvement",
                              "Test Phase", "Car Ban", "Energy Label", "Car Tax",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Road Pricing x Campaign", "Road Pricing x Expert Involvement",
                              "Road Pricing x Test Phase", "Campaign x Test Phase", "Campaign x Expert Involvement",
                              "Test Phase x Expert Involvement", "Road Pricing x Campaign x Test Phase",
                              "Road Pricing x Campaign x Expert Involvement", "Road Pricing x Test Phase x Expert Involvement",
                              "Campaign x Test Phase x Expert Involvement", "Road Pricing x Campaign x Test Phase x Expert Involvement"),
        custom.model.names = c("Two-Way Interaction", "Three Way Interaction", "Four Way Interaction"),
        stars = c(0.01, 0.05, 0.1))

############ Table A5 #############
mCB_choiceI2 <- lmer(choice ~ attrib1_lab *attrib6_lab + attrib1_lab* attrib8_lab + attrib1_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mCB_choiceI3 <- lmer(choice ~ attrib1_lab *attrib6_lab * attrib8_lab +
                       attrib1_lab * attrib6_lab * attrib9_lab +
                       attrib1_lab * attrib8_lab * attrib9_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mCB_choiceI4 <- lmer(choice ~ attrib1_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

htmlreg(list(mCB_choiceI2, mCB_choiceI3, mCB_choiceI4),
        file = "./tables/JPP/A5_PackagingCB_Choice.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Car Ban", "Campaign", "Expert Involvement",
                              "Test Phase", "Energy Label", "Road Pricing", "Car Tax",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Car Ban x Campaign", "Car Ban x Expert Involvement",
                              "Car Ban x Test Phase", "Campaign x Test Phase", "Campaign x Expert Involvement",
                              "Test Phase x Expert Involvement", "Car Ban x Campaign x Test Phase",
                              "Car Ban x Campaign x Expert Involvement", "Car Ban x Test Phase x Expert Involvement",
                              "Campaign x Test Phase x Expert Involvement", "Car Ban x Campaign x Test Phase x Expert Involvement"),
        custom.model.names = c("Two-Way Interaction", "Three Way Interaction", "Four Way Interaction"),
        stars = c(0.01, 0.05, 0.1))

############ Table A6 #############
mRP_choiceI2 <- lmer(choice ~ attrib3_lab *attrib6_lab + attrib3_lab * attrib9_lab + attrib3_lab* attrib8_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mRP_choiceI3 <- lmer(choice ~ attrib3_lab *attrib6_lab * attrib8_lab+ attrib3_lab *attrib6_lab * attrib9_lab + attrib3_lab* attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mRP_choiceI4 <- lmer(choice ~ attrib3_lab *attrib6_lab * attrib8_lab*  attrib9_lab  + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

htmlreg(list(mRP_choiceI2, mRP_choiceI3, mRP_choiceI4),
        file = "./tables/JPP/A6_PackagingRP_Choice.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Road Pricing", "Campaign", "Expert Involvement",
                              "Test Phase", "Car Ban", "Energy Label", "Car Tax",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Road Pricing x Campaign", "Road Pricing x Expert Involvement",
                              "Road Pricing x Test Phase", "Campaign x Test Phase", "Campaign x Expert Involvement",
                              "Test Phase x Expert Involvement", "Road Pricing x Campaign x Test Phase",
                              "Road Pricing x Campaign x Expert Involvement", "Road Pricing x Test Phase x Expert Involvement",
                              "Campaign x Test Phase x Expert Involvement", "Road Pricing x Campaign x Test Phase x Expert Involvement"),
        custom.model.names = c("Two-Way Interaction", "Three Way Interaction", "Four Way Interaction"),
        stars = c(0.01, 0.05, 0.1))

############ Table A7 #############
mCT_choiceI2 <- lmer(choice ~ attrib4_lab *attrib6_lab + attrib4_lab* attrib8_lab + attrib4_lab* attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mCT_choiceI3 <- lmer(choice ~ attrib4_lab *attrib6_lab * attrib8_lab+
                       attrib4_lab *attrib6_lab * attrib9_lab +
                       attrib4_lab *attrib9_lab * attrib8_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mCT_choiceI4 <- lmer(choice ~ attrib4_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

htmlreg(list(mCT_choiceI2, mCT_choiceI3, mCT_choiceI4),
        file = "./tables/JPP/A7_PackagingCT_Choice.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Car Tax", "Campaign", "Expert Involvement",
                              "Test Phase", "Car Ban", "Energy Label", "Road Pricing",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Car Tax x Campaign", "Car Tax x Expert Involvement",
                              "Car Tax x Test Phase", "Campaign x Test Phase", "Campaign x Expert Involvement",
                              "Test Phase x Expert Involvement", "Car Tax x Campaign x Test Phase",
                              "Car Tax x Campaign x Expert Involvement", "Car Tax x Test Phase x Expert Involvement",
                              "Campaign x Test Phase x Expert Involvement", "Car Tax x Campaign x Test Phase x Expert Involvement"),
        custom.model.names = c("Two-Way Interaction", "Three Way Interaction", "Four Way Interaction"),
        stars = c(0.01, 0.05, 0.1))

############ Table A8 #############
mEL_choiceI2 <- lmer(choice ~ attrib2_lab *attrib6_lab + attrib2_lab * attrib8_lab + attrib2_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mEL_choiceI3 <- lmer(choice ~ attrib2_lab *attrib6_lab * attrib8_lab +
                       attrib2_lab *attrib6_lab * attrib9_lab +
                       attrib2_lab *attrib9_lab * attrib8_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)
mEL_choiceI4 <- lmer(choice ~ attrib2_lab *attrib6_lab * attrib8_lab * attrib9_lab + attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab +  attrib7_lab + (1 | id), data = df_conj_out2)

htmlreg(list(mEL_choiceI2, mEL_choiceI3, mEL_choiceI4),
        file = "./tables/JPP/A8_PackagingEL_Choice.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Energy Label", "Campaign", "Expert Involvement",
                              "Test Phase", "Car Ban", "Road Pricing", "Car Tax",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Energy Label x Campaign", "Energy Label x Expert Involvement",
                              "Energy Label x Test Phase", "Campaign x Test Phase", "Campaign x Expert Involvement",
                              "Test Phase x Expert Involvement", "Energy Label x Campaign x Test Phase",
                              "Energy Label x Campaign x Expert Involvement", "Energy Label x Test Phase x Expert Involvement",
                              "Campaign x Test Phase x Expert Involvement", "Energy Label x Campaign x Test Phase x Expert Involvement"),
        custom.model.names = c("Two-Way Interaction", "Three Way Interaction", "Four Way Interaction"),
        stars = c(0.01, 0.05, 0.1))

############ Table A9 #############
m1_sum_choice <- lmer(choice ~ sumpm + sumsm + (1 | id), data = df_conj_out2)
m1_sum_rate <- lmer(rate ~ sumpm + sumsm + (1 | id), data = df_conj_out2)
m2_sum_rate <- lmer(rate ~ sumpm + sumsm + frame + age + gender + edu + usagecar + (1 | id), data = df_conj_out2)
m5_sum_choice <- lmer(choice ~ sumpm + sumsm + sumpm_sq + sumsm_sq + sumpm_sumsm + (1 | id), data = df_conj_out2)
m6_sum_rate <- lmer(rate ~ sumpm + sumsm + sumpm_sq + sumsm_sq + sumpm_sumsm + frame + age + gender + edu + usagecar + (1 | id), data = df_conj_out2)

htmlreg(list(m1_sum_choice, m1_sum_rate, m2_sum_rate, m5_sum_choice, m6_sum_rate), file = "./tables/JPP/A9_RegressionConj2SumOutjoint.html",
        leading.zero = T,
        digits = 2,
        custom.coef.map = list("(Intercept)"="Intercept", 
                               "sumpm"="No. of Primary Meas.", 
                               "sumsm"="No. of Ancillary Meas.", 
                               "frameEmissions"="Emission Frame", 
                               "age"="Age", 
                               "gendermale"="Male", 
                               "eduOther"="Other education", 
                               "eduSecondary Education"="Secondary education", 
                               "eduTertiary Education"="Tertiary education", 
                               "usagecar"="Weekly car usage",
                               "sumpm_sumsm"="Interaction", 
                               "sumpm_sq"="Primary Meas. sq.", 
                               "sumsm_sq"="Ancillary Meas. sq."),
        reorder.coef = c(2:3,11:13,4:10,1),
        custom.model.names = c("Choice: Baseline", "Rate: Baseline", "Rate: Controls", "Choice: Full", "Rate: Full"),
        stars = c(0.01, 0.05, 0.1))

############ Table A10 #############
m_rate <- lmer(rate ~ attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab + attrib6_lab + attrib7_lab + attrib8_lab + attrib9_lab + (1 | id), data = df_conj_out2)
m_choice <- lmer(choice ~ attrib1_lab + attrib2_lab + attrib3_lab + attrib4_lab + attrib5_lab + attrib6_lab + attrib7_lab + attrib8_lab + attrib9_lab + (1 | id), data = df_conj_out2)
m_prim_rate <- lmer(rate ~ attrib1_lab * attrib2_lab * attrib3_lab * attrib4_lab + attrib5_lab + attrib6_lab + attrib7_lab + attrib8_lab + attrib9_lab + (1 | id), data = df_conj_out2)
m_prim_choice <- lmer(choice ~ attrib1_lab * attrib2_lab * attrib3_lab * attrib4_lab + attrib5_lab + attrib6_lab + attrib7_lab + attrib8_lab + attrib9_lab + (1 | id), data = df_conj_out2)

htmlreg(list(m_choice, m_prim_choice,m_rate, m_prim_rate),
        file = "./tables//JPP/A10_PackagingPrimary.html",
        leading.zero = T,
        digits = 2,
        single.row = T,
        custom.coef.names = c(NA, "Car Ban", "Energy Label", "Road Pricing", "Car Tax",
                              "Revenue Use - Street Bulding, Maintenance & Public Transport",
                              "Revenue Use - General Budget", "Public Transportation Price Reduction",
                              "Campaign","Implementation 2030", "Implementation 2035",
                              "Implementation 2040", "Implementation 2045",
                              "Test Phase","Expert Involvement",
                              "Car Ban x Energy Label", "Car Ban x Road Pricing",
                              "Energy Label x Road Pricing", "Car Ban x Car Tax",
                              "Energy Label x Car Tax", "Road Pricing x Car Tax",
                              "Car Ban x Energy Label x Road Pricing", "Car Ban x Energy Label x Car Tax",
                              "Car Ban x Road Pricing x Car Tax", "Energy Label x Road Pricing x Car Tax",
                              "Car Ban x Energy Label x Road Pricing x Car Tax"),
        custom.model.names = c("Baseline Choice", "4w Interaction Baseline", "Baseline Rating", "4w Interaction Rating"),
        stars = c(0.01, 0.05, 0.1))


## END OF SCRIPT